package analysis;


/**
 * @author Mathieu Rousseau
 * @email mathieu.rousseau3@mail.mcgill.ca
 */


import java.io.BufferedReader;
import java.io.BufferedWriter;
import java.io.FileReader;
import java.io.FileWriter;
import java.io.IOException;
import java.text.DecimalFormat;
import java.text.NumberFormat;
import java.util.ArrayList;
import java.util.Collection;
import java.util.Collections;

import utilities.PDButility;


/** ------------------------------------------------------------------------------------------------------------------
 * start of Clustering class
 * ------------------------------------------------------------------------------------------------------------------  */
public class Clustering {

	
	/*  global variables  */
	private final String inputPath = "Clustering/Input/";
	private final String outputPath = "Clustering/Output/";
	private String experimentName = "no_name";
	private int numberOfStructures = 0;
	private int structureFilenameStartIndex = 1;
	private int numberOfPoints = 0;
	private int intraStructureDistanceMatrixSize = 0;
	private int[][] intraConvertIndex;
	private double[][] intraStructureDistanceMatrix;
	private int interStructureDistanceMatrixSize = 0;
	private int[][] interConvertIndex;
	private double[] interStructureDistanceMatrix;
	private ArrayList<Node> hierarchicalTree;
	private NumberFormat formatter = new DecimalFormat("#0.00");
	private PDButility pdb = new PDButility();
	/*  end of global variables  */
	
	
	/*  start of setupIntraConvertIndexMatrix method  */
	private void setupIntraConvertIndexMatrix(){		
		for(int i=0; i<numberOfPoints-1; i++){
			for(int j=i+1; j<numberOfPoints; j++){
				//	find the index value in the optimized 1D array
				int afterSum = ((numberOfPoints-1-i) * (numberOfPoints-i)) / 2;
				intraConvertIndex[i][j] = intraStructureDistanceMatrixSize - afterSum + ((j - i) - 1); 
			}
		}
	}/*  end of setupIntraConvertIndexMatrix method  */
	
	
	/*  start of setupInterConvertIndexMatrix method  */
	private void setupInterConvertIndexMatrix(){		
		for(int i=0; i<numberOfStructures-1; i++){
			for(int j=i+1; j<numberOfStructures; j++){
				//	find the index value in the optimized 1D array
				int afterSum = ((numberOfStructures-1-i) * (numberOfStructures-i)) / 2;
				interConvertIndex[i][j] = interStructureDistanceMatrixSize - afterSum + ((j - i) - 1); 
			}
		}
	}/*  end of setupInterConvertIndexMatrix method  */
	
	
	/*  start of calculateIntraStructureDistance method  */
	private void calculateIntraStructureDistance(int structureIndex, int intraStructureMatrixIndex){
		//  define the structure filename
		String structureFilename = "default.pdb";
		if(structureIndex <= 9){
			structureFilename = "structure_000" + structureIndex + "_" + experimentName + ".pdb";
		}
		else if(structureIndex <= 99){
			structureFilename = "structure_00" + structureIndex + "_" + experimentName + ".pdb";
		}
		else if(structureIndex <= 999){
			structureFilename = "structure_0" + structureIndex + "_" + experimentName + ".pdb";
		}
		else{
			structureFilename = "structure_" + structureIndex + "_" + experimentName + ".pdb";
		}
		
		
		//  get the structure from file
		double[][] structure = pdb.getStructure(inputPath + structureFilename);
		
		
		//  calculate and store the intra structure distances
		for(int i=0; i<numberOfPoints-1; i++){
			for(int j=i+1; j<numberOfPoints; j++){
				double term1 = structure[i][0] - structure[j][0];
				term1 *= term1;
				double term2 = structure[i][1] - structure[j][1];
				term2 *= term2;
				double term3 = structure[i][2] - structure[j][2];
				term3 *= term3;
				double distance = Math.sqrt(term1 + term2 + term3);
				
				intraStructureDistanceMatrix[intraStructureMatrixIndex][intraConvertIndex[i][j]] = distance;
			}
		}
		
		
		//  make sure that the memory used by the structure is freed
		structure = null;
	}/*  end of calculateIntraStructureDistance method  */
	
	
	/*  start of calculateInterStructureDistance method  */
	private double calculateInterStructureDistance(int structureIndex1, int structureIndex2){
		double distanceBetweenStructures = 0.0;
		for(int i=0; i<numberOfPoints-1; i++){
			for(int j=i+1; j<numberOfPoints; j++){
				double difference = intraStructureDistanceMatrix[structureIndex1][intraConvertIndex[i][j]] - intraStructureDistanceMatrix[structureIndex2][intraConvertIndex[i][j]];
				double maxDistance = Math.max(intraStructureDistanceMatrix[structureIndex1][intraConvertIndex[i][j]], intraStructureDistanceMatrix[structureIndex2][intraConvertIndex[i][j]]);				
				distanceBetweenStructures += Math.pow(difference / maxDistance, 2);
			}
		}
		return distanceBetweenStructures;
	}/*  end of calculateInterStructureDistance method  */
	
	
	/*  start of printIntraConvertIndexMatrix method  */
	private void printIntraConvertIndexMatrix(){
		for(int i=0; i<numberOfPoints; i++){
			for(int j=0; j<numberOfPoints; j++){
				if(j<=i){
					System.out.print("\t");
				}
				else{
					System.out.print(intraConvertIndex[i][j] + "\t");
				}
			}
			System.out.println();
		}
	}/*  end of printIntraConvertIndexMatrix method  */
	

	/*  start of printInterConvertIndexMatrix method  */
	private void printInterConvertIndexMatrix(){
		for(int i=0; i<numberOfStructures; i++){
			for(int j=0; j<numberOfStructures; j++){
				if(j<=i){
					System.out.print("\t");
				}
				else{
					System.out.print(interConvertIndex[i][j] + "\t");
				}
			}
			System.out.println();
		}
	}/*  end of printInterConvertIndexMatrix method  */
	
	
	/*  start of printIntraDistanceMatrix method  */
	private void printIntraDistanceMatrix(int structureIndex){
		for(int i=0; i<numberOfPoints; i++){
			for(int j=0; j<numberOfPoints; j++){
				if(j<=i){
					System.out.print("\t");
				}
				else{
					System.out.print(formatter.format(intraStructureDistanceMatrix[structureIndex][intraConvertIndex[i][j]]) + "\t");
				}				
			}
			System.out.println();
		}
	}/*  end of printIntraDistanceMatrix method  */
	
	
	/*  start of printInterDistanceMatrix method  */
	private void printInterDistanceMatrix(){
		for(int i=0; i<numberOfStructures; i++){
			for(int j=0; j<numberOfStructures; j++){
				if(j<=i){
					System.out.print("\t");
				}
				else{
					System.out.print(formatter.format(interStructureDistanceMatrix[interConvertIndex[i][j]]) + "\t");
				}				
			}
			System.out.println();
		}
	}/*  end of printInterDistanceMatrix method  */
	
	
	/*  start of outputInterDistanceMatrixToFile method  */
	private void outputInterDistanceMatrixToFile(){
		try{						
			BufferedWriter bw = new BufferedWriter(new FileWriter(outputPath + "interStructureDistanceMatrix.txt"));
			
			
			//  write the column headers
			for(int i=0; i<numberOfStructures; i++){
				if(i<numberOfStructures-1){
					bw.write((i+1) + "\t");
				}
				else{
					bw.write(String.valueOf(i+1));
				}
			}
			bw.flush();
			bw.newLine();
			
			
			//  write the distance values
			for(int i=0; i<numberOfStructures; i++){
				//  write the row label
				bw.write((i+1) + "\t");
				
				for(int j=0; j<numberOfStructures; j++){
					if((i == numberOfStructures-1) && (j == numberOfStructures-1)){
						bw.write("0.00");		
						bw.flush();
					}
					else if(i == j){
						bw.write("0.00\t");
					}
					else if(j<=i){
						bw.write(formatter.format(interStructureDistanceMatrix[interConvertIndex[j][i]]) + "\t");
						bw.flush();
					}
					else if(j<numberOfStructures-1){
						bw.write(formatter.format(interStructureDistanceMatrix[interConvertIndex[i][j]]) + "\t");
						bw.flush();
					}
					else{
						bw.write(formatter.format(interStructureDistanceMatrix[interConvertIndex[i][j]]));
						bw.flush();
					}
				}
				bw.newLine();
				bw.flush();
			}			
			bw.flush();
			bw.close();
		}catch (IOException e) {
			e.printStackTrace();
		}//end of try-catch block
		System.out.println("Output file: interStructureDistanceMatrix.txt");
	}/*  end of outputInterDistanceMatrixToFile method  */
	
	
	/*  start of clusterStructures method  */
	private void clusterStructures(ArrayList<Node> tree){
		
		//  add all of the structures to the list of nodes as leaf nodes
		for(int i=0; i<numberOfStructures; i++){
			ArrayList<Integer> tempList = new ArrayList<Integer>();
			tempList.add(new Integer(i));
			tempList.trimToSize();
			Node leafStructure = new Node(tempList);
			tree.add(leafStructure);
		}
		tree.trimToSize();
		
		
		//  cluster the nodes using hierarchical clustering and average-linkage metric
		while(tree.size() > 1){			
			groupClosestNodes(tree);
		}
	}/*  end of clusterStructures method  */
	
	
	/*  start of groupClosestNodes method  */
	private void groupClosestNodes(ArrayList<Node> tree){
		double minDistance = Double.MAX_VALUE;
		Node[] closestNodes = {null, null};
				
		for(Node node1: tree){
			for(Node node2: tree){
				if(node1 != node2){
				
					//  using average linkage metric
					double averageDistance = 0.0;
					
					for(Integer structureIndexNode1: node1.getListOfStructures()){
						for(Integer structureIndexNode2: node2.getListOfStructures()){
							averageDistance += interStructureDistanceMatrix[interConvertIndex[structureIndexNode1.intValue()][structureIndexNode2.intValue()]];
						}					
					}
					averageDistance = averageDistance / (node1.getNumberOfStructures() * node2.getNumberOfStructures());
					if(averageDistance < minDistance){
						minDistance = averageDistance;
						closestNodes[0] = node1;
						closestNodes[1] = node2;
					}
				}
			}
		}
		
		
		//  create the new node and set it as the parent node of the two closest nodes
		Node newNode = new Node(closestNodes[0].getListOfStructures(), closestNodes[1].getListOfStructures(), null, closestNodes[0], closestNodes[1], minDistance - getTreeHeight(closestNodes[0]), minDistance - getTreeHeight(closestNodes[1]));
		closestNodes[0].setParent(newNode);
		closestNodes[1].setParent(newNode);
		
				
		//  remove the node with the biggest index first to avoid list-shift problems
		for(Node nodeElement: tree){
			if(closestNodes[0] == nodeElement){
				tree.remove(closestNodes[1]);
				tree.remove(closestNodes[0]);
				break;
			}
			else if(closestNodes[1] == nodeElement){
				tree.remove(closestNodes[0]);
				tree.remove(closestNodes[1]);
				break;
			}
		}
		
		
		//  add the new node to the tree
		tree.add(newNode);
		
		
		//  trim the tree
		tree.trimToSize();
	}/*  end of groupClosestNodes method  */
	
	
	/*  start of getTreeHeight method  */
	private double getTreeHeight(Node node){
		//  the height of the left subtree is the same as the height of the right subtree
		if(node.getLeftChild() == null){
			return 0;
		}
		else{
			return node.getLeftBranchLength() + getTreeHeight(node.getLeftChild());
		}
	}/*  end of getTreeHeight method  */
	
	
	/*  start of printTree method  */
	private void printTree(){
		for(Node nodeElement: hierarchicalTree){
			System.out.print("Number of structures = " + nodeElement.getNumberOfStructures() + ":\t");
			for(Integer structureIndex: nodeElement.getListOfStructures()){
				System.out.print(structureIndex.intValue() + "\t");
			}
			System.out.println();
		}
	}/*  end of printTree method  */
	
	
	/*  start of outputTreeToFile method  */
	private void outputTreeToFile(Node root){
		if((root.getLeftChild() == null) && (root.getRightChild() == null)){
			System.out.println("Error: the tree only has a root node. There are no children nodes.\n" +
					"Cannot output the tree to file.");
		}
		else{
			try{						
				BufferedWriter bw = new BufferedWriter(new FileWriter(outputPath + "hierarchicalTree.dendro"));
			
				String NewickFormatOutputString = new String();
				
				NewickFormatOutputString = recursiveTreeTraversal(root);
				NewickFormatOutputString += ";";
				
//				System.out.println(NewickFormatOutputString);
				
				bw.write(NewickFormatOutputString);
				bw.newLine();
				
				bw.flush();
				bw.close();
			}catch (IOException e) {
				e.printStackTrace();
			}//end of try-catch block
			System.out.println("Output file: hierarchicalTree.dendro");
		}
	}/*  end of outputTreeToFile method  */
	
	
	/*  start of recursiveTreeTraversal method  */
	private String recursiveTreeTraversal(Node node){
		String returnString = new String();
		if(node.getNumberOfStructures() == 1){
			return String.valueOf((node.getListOfStructures().get(0).intValue()+1));			//  start the structure indices from 1 when printing
		}
		else{
			returnString += "(";
			returnString += recursiveTreeTraversal(node.getLeftChild());
			returnString += ":";
			returnString += formatter.format(node.getLeftBranchLength());
			returnString += ",";
			returnString += recursiveTreeTraversal(node.getRightChild());
			returnString += ":";
			returnString += formatter.format(node.getRightBranchLength());
			returnString += ")";
		}
		return returnString;
	}/*  end of recursiveTreeTraversal method  */
	
	
	/*  start of run method  */
	public void run(String[] args) {
		if(args.length == 2){
			experimentName = args[0];
			numberOfStructures = Integer.parseInt(args[1]);
		}
		else if(args.length == 3){
			experimentName = args[0];
			numberOfStructures = Integer.parseInt(args[1]);
			structureFilenameStartIndex = Integer.parseInt(args[2]);
		}
		else{
//			System.out.println("Error:\tplease use one of the following patterns (without the quotes)\n" +
//					"\t\"java ClusteringMain experimentName numberOfStructures\"\n" +
//					"\t\"java ClusteringMain experimentName numberOfStructures structureFilenameStartIndex\"\n" +					
//				"Notes:\tthe default experiment name is " + experimentName + "\n" +
//					"\tthe default number of structures to cluster is " + numberOfStructures + "\n" +
//					"\tthe default structureFilenameStartIndex is " + structureFilenameStartIndex + "\n" +
//				"System.exit(0)");
//			System.exit(0);
		}
		
		
		
		//  set the name of the experiment which is present in all of the structure file names
		//  format for the structure filenames: "structure_####_experimentName.pdb"
		experimentName = "MCMC5C";
		
		
		//  set the number of structures that will be clustered
		numberOfStructures = 400;
		
		
		//  set the start index of the structures that will be clustered
		structureFilenameStartIndex = 1;
		
		
		//  find out the number of points in the structures
		numberOfPoints = pdb.getNumberOfPoints(inputPath + "structure_0001_" + experimentName + ".pdb");
		
		
		//  setup the intra convert index matrix for memory usage optimization
		intraStructureDistanceMatrixSize = (numberOfPoints * (numberOfPoints-1)) / 2;
		intraConvertIndex = new int[numberOfPoints][numberOfPoints];
		setupIntraConvertIndexMatrix();
//		printIntraConvertIndexMatrix();
		
		
		//  initialize the matrix to store the intra distances for all the structures
		intraStructureDistanceMatrix = new double[numberOfStructures][intraStructureDistanceMatrixSize];
		
		
		//  convert all of the structures to intra distance matrices
		for(int i=structureFilenameStartIndex; i<(structureFilenameStartIndex+numberOfStructures); i++){
			calculateIntraStructureDistance(i, i-structureFilenameStartIndex);
		}
//		printIntraDistanceMatrix(0);
//		for(int i=0; i<numberOfStructures; i++){
//			printIntraDistanceMatrix(i);
//			System.out.println("\n");
//		}
		
		
		//  setup the inter convert index matrix for memory usage optimization
		interStructureDistanceMatrixSize = (numberOfStructures * (numberOfStructures - 1)) / 2;
		interConvertIndex = new int[numberOfStructures][numberOfStructures];
		setupInterConvertIndexMatrix();
//		printInterConvertIndexMatrix();
		
		
		//  initialize the matrix to store the inter distances for all the structures
		interStructureDistanceMatrix = new double[interStructureDistanceMatrixSize];

		
		//  calculate the distance between all of the structures
		for(int i=0; i<numberOfStructures-1; i++){
			for(int j=i+1; j<numberOfStructures; j++){
				interStructureDistanceMatrix[interConvertIndex[i][j]] = calculateInterStructureDistance(i, j);
			}
		}
//		printInterDistanceMatrix();
		outputInterDistanceMatrixToFile();
				
		
		//  cluster the structures hierarchically
		hierarchicalTree = new ArrayList<Node>();
		clusterStructures(hierarchicalTree);
//		printTree();
		outputTreeToFile(hierarchicalTree.get(0));				//  the final list only contains the root node of the tree
		
		
	}/*  end of run method  */

	
	
	
	/** ------------------------------------------------------------------------------------------------------------------
	 * start of Node class
	 * ------------------------------------------------------------------------------------------------------------------  */
	private class Node{
		
		
		/*  node fields  */
		private ArrayList<Integer> listOfStructures;
		private double leftBranchLength;
		private double rightBranchLength;
		private Node parent;
		private Node leftChild;
		private Node rightChild; 
		/*  end of node fields  */
		
		
		/*  start of first constructor method  */
		Node(){
			this(null, null, null, null, null, 0, 0); 
		}/*  end of first constructor method  */
		
		
		/*  start of second constructor method  */
		Node(ArrayList<Integer> structureIndexList){
			this(structureIndexList, null, null, null, null, 0, 0);
		}/*  end of second constructor method  */
		
		
		/*  start of third constructor method  */
		Node(ArrayList<Integer> structureList1, ArrayList<Integer> structureList2, Node parent, Node leftChild, Node rightChild, double leftBranchLength, double rightBranchLength){
			//  add all of the structure indices to the list
			listOfStructures = new ArrayList<Integer>();
			if(structureList1 != null){
				for(Integer structureIndex:structureList1){
					Integer structureIndexDeepCopy = new Integer(structureIndex.intValue());
					listOfStructures.add(structureIndexDeepCopy);
				}
			}
			if(structureList2 != null){
				for(Integer structureIndex:structureList2){
					Integer structureIndexDeepCopy = new Integer(structureIndex.intValue());
					listOfStructures.add(structureIndexDeepCopy);
				}
			}
			listOfStructures.trimToSize();
			
			
			//  set the parent and child nodes
			this.parent = parent;
			this.leftChild = leftChild;
			this.rightChild = rightChild;
			
			
			//  set the branch lengths to the children
			this.leftBranchLength = leftBranchLength;
			this.rightBranchLength = rightBranchLength;
			//  check to make sure that the branch lengths are not negative
			if(this.leftBranchLength < 0){
				System.out.println("Error: trying to create a node with a negative (" + formatter.format(this.leftBranchLength) + ") left branch length.\n" +
						"Setting the left branch length to zero (0)");
				this.leftBranchLength = 0;
			}
			if(this.rightBranchLength < 0){
				System.out.println("Error: trying to create a node with a negative (" + formatter.format(this.rightBranchLength) + ") right branch length.\n" +
						"Setting the left branch length to zero (0)");
				this.rightBranchLength = 0;
			}
		}/*  end of third constructor method  */
		
		
		/*  start of getNumberOfStructures method  */
		public int getNumberOfStructures(){
			listOfStructures.trimToSize();
			return listOfStructures.size();
		}/*  end of getNumberOfStructures method  */


		/*  start of getListOfStructures method  */
		public ArrayList<Integer> getListOfStructures() {
			return listOfStructures;
		}/*  end of getListOfStructures method  */


		/*  start of setListOfStructures method  */
		public void setListOfStructures(ArrayList<Integer> listOfStructures) {
			this.listOfStructures = listOfStructures;
		}/*  end of setListOfStructures method  */
		
		
		/*  start of setListOfStructures method  */
		public void setListOfStructures(ArrayList<Integer> structureList1, ArrayList<Integer> structureList2) {
			//  add all of the structure indices to the list
			this.listOfStructures = new ArrayList<Integer>();
			for(Integer structureIndex:structureList1){
				Integer structureIndexDeepCopy = new Integer(structureIndex.intValue());
				this.listOfStructures.add(structureIndexDeepCopy);
			}
			for(Integer structureIndex:structureList2){
				Integer structureIndexDeepCopy = new Integer(structureIndex.intValue());
				this.listOfStructures.add(structureIndexDeepCopy);
			}
			this.listOfStructures.trimToSize();			
		}/*  end of setListOfStructures method  */


		/*  start of getLeftBranchLength method  */
		public double getLeftBranchLength() {
			return leftBranchLength;
		}/*  end of getLeftBranchLength method  */


		/*  start of setLeftBranchLength method  */
		public void setLeftBranchLength(double leftBranchLength) {
			this.leftBranchLength = leftBranchLength;
		}/*  end of setLeftBranchLength method  */


		/*  start of getRightBranchLength method  */
		public double getRightBranchLength() {
			return rightBranchLength;
		}/*  end of getRightBranchLength method  */


		/*  start of setRightBranchLength method  */
		public void setRightBranchLength(double rightBranchLength) {
			this.rightBranchLength = rightBranchLength;
		}/*  end of setRightBranchLength method  */


		/*  start of getParent method  */
		public Node getParent() {
			return parent;
		}/*  end of getParent method  */


		/*  start of setParent method  */
		public void setParent(Node parent) {
			this.parent = parent;
		}/*  end of setParent method  */


		/*  start of getLeftChild method  */
		public Node getLeftChild() {
			return leftChild;
		}/*  end of getLeftChild method  */


		/*  start of setLeftChild method  */
		public void setLeftChild(Node leftChild) {
			this.leftChild = leftChild;
		}/*  end of setLeftChild method  */


		/*  start of getRightChild method  */
		public Node getRightChild() {
			return rightChild;
		}/*  end of getRightChild method  */


		/*  start of setRightChild method  */
		public void setRightChild(Node rightChild) {
			this.rightChild = rightChild;
		}/*  end of setRightChild method  */
		
		
		/*  start of toString override method  */
		public String toString(){
			String nodeString = "Number of structures = " + this.getNumberOfStructures() + "\n";
			if(this.getNumberOfStructures() == 1){
				nodeString += "LeafNode:\t" + this.getListOfStructures().get(0).intValue();
			}
			else{
				Node leftChild = this.getLeftChild();
				nodeString += "LeftChild: (n=" + leftChild.getNumberOfStructures() + " | d=" + formatter.format(this.getLeftBranchLength()) + ")\t";				
				for(Integer structureElement: leftChild.getListOfStructures()){
					nodeString += structureElement.intValue() + "\t";
				}
				nodeString += "\n";
				
				Node rightChild = this.getRightChild();
				nodeString += "RightChild: (n=" + rightChild.getNumberOfStructures() + " | d=" + formatter.format(this.getRightBranchLength()) + ")\t";				
				for(Integer structureElement: rightChild.getListOfStructures()){
					nodeString += structureElement.intValue() + "\t";
				}
			}
			return nodeString;
		}/*  end of toString method  */
		
		
	}/*  end of Node class  */

	
}/*  end of Clustering class  */
